Microstructure identification via detrended fluctuation analysis of ultrasound signals 
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We describe an algorithm for simulating ultrasound propagation in random one-dimensional media, mim- 
icking different microstructures by choosing physical properties such as domain sizes and mass densities from 
probability distributions. By combining a detrended fluctuation analysis (DFA) of the simulated ultrasound sig- 
nals with tools from the pattern-recognition literature, we build a Gaussian classifier which is able to associate 
each ultrasound signal with its corresponding microstructure with a very high success rate. Furthermore, we also 
show that DFA data can be used to train a multilayer perceptron which estimates numerical values of physical 
properties associated with distinct microstructures. 



I. INTRODUCTION 

Much attention has been given to the problem of wave prop- 
agation in random media by the condensed-matter physics 
community, especially in the context of Anderson localization 
and its analogs 1 1 -5 1. A hallmark of these phenomena is the 
fact that randomness induces wave attenuation by energy con- 
finement, even in cases where dissipation can be neglected. 
The interplay of energy confinement and randomness gives 
rise to noisy but correlated signals, for instance, in the case of 
acoustic pulses propagating in inhomogeneous media. 

It has long been known that the correlations in a time series 
hide relevant information on the dynamics which generated 
it. In a pioneering paper, Hurst [6] introduced the rescaled- 
range analysis of a time series, which measures the power- 
law growth of properly rescaled flucuations in the series as 
one looks at larger and larger time scales x. The associated 
Hurst exponent H governing the growth of such fluctuations 
is able to gauge memory effects on the underlying dynamics, 
offering insight into its character. The presence of long-term 
memory leads to a value of the exponent H which deviates 
from the uncorrected random-walk value H = 1/2, persistent 
(antipersistent) behavior of the time series yielding H > 1/2 
(H < 1/2). Additionally, a crossover at a time scale T x be- 
tween two regimes characterized by different Hurst exponents 
may reveal the existence of competing ingredients in the dy- 
namics, and in principle provides a signature of the associated 
system. This forms the base of methods designed to charac- 
terize such systems [7|, if one is able to obtain reliable esti- 
mates of the Hurst exponents. This is in general a difficult 
task, as even local trends superposed on the noisy signal may 
affect the rescaled-range analysis, obscuring the value of H. 
A related exponent, a, defined through detrended flucuation 
analysis (DFA) [8|, can be used instead. 

Actually, the characteristics of exponents and crossovers 
observed in the DFA curves associated with various types of 



data series have been extensively used to distinguish between 
the systems producing such series. Examples include cod- 
ing versus noncoding DNA regions [8 9|, healthy versus dis- 
eased subjects as regards cardiac IflOl . neurological IfTTl [T2l 
and respiratory function fl3l . and ocean versus land regions 
as regards temperature variations 1141 . However, there are 
many instances in which these characteristics are not clearcut 
and the DFA curves exhibit a more complex dependence on 
the time scale. In such cases, it has been shown that pattern 
recognition tools |[T5l can help the identification of relevant 
features, greatly improving the success of classification tasks 

umbo. 

In this work, we investigate the possibility of extracting 
information on the nature of inhomogeneities by analyzing 
fluctuations in time series associated with ultrasound propaga- 
tion in random media. A hint that this possibility is real was 
provided by the fact that the crossover features of DFA and 
Hurst analysis curves from backscattered ultrasound signals 
revealed signatures of the microstructure of cast-iron sam- 
ples J7J. Here, in order to perform a systematic study, we 
resort to simulating the propagation of ultrasound pulses in 
one-dimensional media with distinct microstructures, defined 
by probability distributions of physical properties such as do- 
main size, density and sound velocity. Although this choice of 
geometry cannot allow for the full phenomenology of sound 
propagation (such as mode conversion from transverse to lon- 
gitudinal sound waves), it makes it possible to generate large 
quantities of simulated data, which are important in order to 
assess the generalizability of our reported results. Moreover, it 
approximately describes normal propagation of sound waves 
in layered media. 

The paper is organized as follows. In Sec. |ll]we present the 
artificial microstructures we produced, and present a sketch of 
the simulation technique; a detailed description is relegated to 
Appendix [A] In Sec. Ill we describe the method of detrended 
fluctuation analysis and its results when applied to our simu- 
lated signals. Then, in Sec. IV we report on an automated 
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classifier which is able to associate, with a very high success 
rate, the DFA curves with the corresponding microstructure. 
Furthermore, we show in Sec. [VJthat the DFA curves can be 
used to train a neural network which predicts numerical val- 



transductor 



-W- 



Microstructure Average domain size (m) Average density (kg/nr 3 ) 







Pi 




p 2 






C\ 









X 


1 


X 


2 X 



Figure 1 : Sketch of the geometry used in simulating ultrasound prop- 
agation in an inhomogeneous medium. See main text for labels. 



ues of physical properties associated with diferent microstruc- 



tures. We close the paper presenting a summary in Sec. VI 



II. SIMULATING ULTRASOUND PROPAGATION 

We are interested in studying ultrasound propagation along 
a one-dimensional medium of width W, with a pulse gener- 
ated in a transductor located in one end of the system. Since 
the medium consists of many different domains, with possibly 
different physical properties (density and sound velocity), in 
general the pulse will be scattered as it propagates towards the 
opposite end, where it will be reflected. Information about the 
microstructure is in principle hidden in the scattered signal, 
which is registred in the transductor as it arrives. 

The domains are labelled by an index j, so that domain j 
extends between xj and xj+i, and is characterized by its den- 
sity pj and its sound velocity Cj, (See Figure fTI) For the 
one-dimensional geometry employed here, the solution of the 
wave equation can be carried out semi-analytically, as detailed 
in Appendix A. For a given choice of physical parameters 
{pj,Cj} of the various domains, the displacement field inside 
the medium, as a function of position x and time t, can be 
written as 



4> (x, t) = Y^tykX-k (x) cos (CO**) , 



(1) 



where k labels the different eigenfrequencies <a k , the function 
X k (x) is explicitly given by 

X k (x) = A jk cos ((O k x/cj)+Bj k sin((O k x/cj) , 

with j such that x> < x < and the coefficients Aj k and 
Bj k are determined from boundary conditions at the interfaces 
separating contiguous domains, while the weights § k are de- 
rived from the initial condition 4>(;t:,0). Here we choose an 
initial pulse contained entirely inside the transductor, with a 
form given by 



4>(x,0) = 4>oe 



sin 



2%f (x — JCtrans) 



c trans 



for x inside the transductor, and <I> (x, 0) = otherwise, where 
<t>o is an amplitude, / is the reference frequency of the trans- 
ductor, Ctrans is the sound velocity inside the transductor, jttrans 
is the position of the left end of the transductor, and y is a 
"damping" factor, introduced so as to make the pulse resem- 
ble those produced by a real transductor. In this work, we 
use / = 10 MHz, c trans = 5800 m/s, and y= 129.31 m" 1 , for a 
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Table I: Average values of physical properties for the 16 different 
microstructures used in this work. The sound velocity is fixed at 
5900 m/s in all cases. 
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Figure 2: Typical simulated ultrasound signals from microstructures 
1 (a) to 4 (d). 



transductor of length 2.32cm, in which 4 wavelengths of the 
pulse can fit. The density inside the transductor is chosen as 
ptrans = 2600kg/m 3 (about the density of quartz). We take 
into account in Eq. ([T]) all eigenfrequencies co^ smaller than 
tt>max = 16 x 2%f, corresponding to 16 times the reference an- 
gular frequency of the transductor. We checked that halving 
the value of C0 max has no relevant effect on the results we re- 
port below. 

The ultrasound signals we keep correspond to time series 
{m, } of the displacement captured at the right end of the trans- 
ductor, with a sampling rate of 50 MHz. Each time series con- 
tains 2048 points, corresponding to about 4 x 10~ 5 seconds. 
As expected, results comparable to the ones described below 
are obtained if one makes use of time series defined by the 
sound pressure at the right end of the transductor. 

We work here with 16 different choices of microstructure, 
combining 4 different average domain sizes and 4 different 
average densities, with characteristics detailed in Table [I] The 
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actual size and density of a domain are chosen from Gaussian 
probability distributions with a standard deviation of 10% the 
average value for the size, and of 2.5% the average value for 
the density. For each of the 16 microstructures, we obtained 
signals from 100 different disorder realizations, with a total of 
1600 signals. Our aim is to identify the microstructure based 
on the analysis of the ultrasound signal. Since we keep the 
average system size fixed at W = 2cm, we assume the same 
sound velocity (5900 m/s) for all microstructures, so that the 
time intervals between signal peaks do not trivially reveal the 
microstructure type. Typical signals are shown in Fig. [2] No- 
tice that fluctuations in the signals increase from microstruc- 
ture 1 to microstructure 2, and decrease for microstructures 
3 to 4. This nonmonotonic behavior of the fluctuations as a 
function of average domain size is due to the fact that the av- 
erage time needed for the wave to propagate through a domain 



in microstructure 2 is about 1/ f 
scattering. 



10 s, thus maximizing 



III. DETRENDEND FLUCTUATION ANALYSIS 

The detrended fluctuation analysis (DFA), introduced by 
Peng et al. [8 1, calculates the detrended fluctuations in a time 
series as a function of a time-window size x. The detrending 
comes from fitting the integrated time series inside each time 
window to a polynomial, and calculating the average variance 
of the residuals. Explicitly, the method works as follows. A 
time series {m,} of length M is initially integrated, yielding a 
new time series yj, 



yj 



the average (u) being taken over all points, 



M 



x ' M 



For each time window 4 of size T, the points inside 4 are fitted 
by a polynomial of degree I (which we take in this work to be 
1 = 1, i.e. a straight line), yielding a local trend y correspond- 
ing to the ordinate of the fit. The variance of the residuals is 
calculated as 



i 



1 iei k 



~ \2 

■yj) > 



and fk (x) is averaged over all intervals to yield the detrended 
fluctuation F (x), 



F(x) 



M 



M — x + 1 being the number of time windows of size x in a 
time series with M points. As defined here, X is the (integer) 
number of points inside a time window, the time increment 
between consecutive points corresponding to the inverse sam- 
pling rate, 2 x 10~ 8 s. 
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Figure 3: Average DFA curves for 8 different microstructures as a 
function of the the time-window size x (offset vertically for clarity). 
The dashed lines have slope 1 /2. 



Notice that here, besides using overlapping time windows, 
we also calculate the variance of the residuals inside each win- 
dow, in a similar spirit to what is done for the detrended cross- 
correlation analysis of Ref. EDI . This approach is slightly 
distinct from the original scheme of Ref. [8], where nonover- 
lapping time windows are employed, and the variance is cal- 
culated for the whole time series. When applied to fractional 
Brownian motion ETTl charaterized by a Hurst exponent H, 
both approaches yield the same exponent within numerical er- 
rors. Interestingly, the performance of the classifier described 



in Sec. IV however, is vastly improved by our approach. 

When applied to a time series generated by a process gov- 
erned by a single dynamics, as for instance in fractional Brow- 
nian motion ET1 . DFA yields a function F (x) following a 
power-law behavior, 

F(x)~Cx a , 

in which C is a constant and a is an exponent which is related 
to the Hurst exponent H, measuring memory effects on the 
dynamics. If persistent (antipersistent) behavior of the time 
series is favored, a is larger (smaller) than 1 /2. 

As shown in Fig. [3] for a subset of 8 microstructures, 
the curves of F (x) calculated from the simulated ultrasonic 
signals do not conform in general to a power-law behavior, 
so that the exponent a is ill-defined, except perhaps for mi- 
crostructures 4 and 8, characterized by the largest average do- 
main sizes, which yield exponents approaching the uncorre- 
lated random-walk value a = 1/2. This same value can be 
approximately identified for the other microstructures if the 
analysis is restricted to a range of time-window sizes x such 
that log 10 x > 2.5, which correspond to time scales greater 
than 6.3 x 10~ 6 s, compatible with the time 6.78 x 10~ 6 s 
needed for the pulse to travel across the medium and return 
to the transductor. At shorter time scales, scattering of the 
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Figure 4: Projection of the DFA vectors on the plane defined by 
the first two principal components of the data, corresponding to the 
eigenvectors of the covariance matrix associated with the two largest 
eigenvalues. Inset: vector components of the eigenvectors for each 
of the 37 directions. 

waves at the interfaces between domains introduces large in- 
terference effects leading to the antipersistent behavior re- 
vealed by the F (x) curves. Such effects, as expected, are 
stronger for microstructures 1, 2, 5, and 6, characterized by 
smaller average domain sizes. Even shorter time- window 
sizes, log 10 T < 0.7, probe time scales inferior to the inverse 
frequency of the pulse, 1/f = 10 _7 s, and, as expected, point 
to a locally persistent behavior of the time series. 

Instead of attempting to manually correlate the signals with 
the various microstructures based on a manual identification 
of the various aspects of the curves, we resort to pattern recog- 
nition tools [15]. To this end, we define for each signal i a 
DFA vector x, whose components correspond to the values of 
the function F (x) at a fixed set {xj} of window sizes. Here 
we build {xj} from the integer part of the integer powers of 
2 1 / 4 , from 4 to 2048, comprising 37 different values of x. 

A visualization of the DFA vectors is hindered by their 
high number of components, n = 37. However, a principal- 
component analysis 1031 can reveal the directions along which 
the data for all 1600 vectors is most scattered. This is done by 
projecting each vector along the principal components, corre- 
sponding to the eigenvectors of the covariance matrix 

£ = ^I>-/z)( Xi -/z) 1 \ 

in which the summation runs over all N = 1600 DFA (column) 
vectors x, and /j, is the average vector, 




The eigenvectors are arranged in decreasing order of their re- 
spective eigenvalues. The first principal component thus cor- 
responds to the direction accounting for the largest variation in 



the data, the second principal component to the second largest 
variation, etc. Figure |4] shows projections of all DFA vec- 
tors on the plane defined by the first two principal compo- 
nents, revealing at the same time a clustering of the data for 
each microstructure and a considerable superposition of the 
data for microstructures which differ only by average density 
(microstructures 1, 5, 9, and 13, for instance). Although this 
superposition is in part an artefact of the two-dimensional pro- 
jection, it is not satisfactorily eliminated when the other direc- 
tions are taken into account. Accordingly, attempts at associ- 
ating a vector x,- to the microstructure whose average vector is 
closer to x, lead to an error rate of about 42%. However, as 
discussed in the next Section, a more sophisticated approach 
considerably improves the classification performance. 

Notice that, as shown in the inset of Fig. |4] special features 
of the first two eigenvectors of the covariance matrix are asso- 
ciated with directions 26 and above, along which the compo- 
nents of the first (second) eigenvector have typically smaller 
(larger) absolute values than along other directions. Moreover, 
the components of the second eigenvector along directions 1 
to 5 also have locally larger absolute values. In fact, direc- 
tion 5 is connected with the inverse frequency 1//, and di- 
rections above 26 are associated with time-window sizes such 
that log 10 T > 2.5, again pointing to the special role played 
by these time scales in differentiating the various microstruc- 
tures. 



IV. GAUSSIAN DISCRIMINANTS 

We first want to check whether it is possible to build an ef- 
ficient automated classifier which is able to assign a signal to 
one of the possible microstructures, based on the correspond- 
ing DFA vector. Attempts at assigning a DFA vector x to the 
microstructure (or class, in the language of pattern recogni- 
tion) whose average vector lies closer to x lead to many mis- 
classifications due to the fact that the average vectors of simi- 
lar microstructures (such as 1 and 5) are close to one another. 
As a classification based solely on distance disregards addi- 
tional information provided by the probability distributions of 
the DFA vectors obtained from each microstructure, here we 
follow an approach to discrimination based on estimates of 
those distributions. 

Our task is to estimate the probability P (coy |x) that a given 
vector x belongs to class C0y, j G { 1 , 2, . . . , C}. (In our case, of 
course, C = 16.) From Bayes' theorem, this probability can 
be written as 

where P(x|cOy) is the probability that a sample belonging to 
class C0y produces a vector x, P (coy) is the prior probability of 
class CO, occurring, and P (x) is the prior probability of vector 
x occurring. Once P((£>j\x) is known for all classes C0y, we 
assign vector x to class C0y if 

P((Oj\x) >P(& k \x) , for all k ^ j. 
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Since P (x) is class-independent, and thus irrelevant to the de- 
cision process, the problem of calculating P(cO/|x) reduces to 
estimating P(x|co,-) and P (CO,-). 

Among the various possibilities for the estimation of 
P (x|(0/), we choose to work with normal-based quadratic dis- 
criminant functions fl5l . derived as follows. We assume that 
P(x|cO;) has the Gaussian form 



P(x\(Oj) = 



1 



(2jc)2|Ej 



exp 



--(x-^fST^X-/*;) 



where n is the number of components of x, while fij and S 7 
are the average vector and the covariance matrix of class CO;. 
By selecting a subset of the available vectors to form a training 
set {x,}, unbiased maximum-likelihood estimates of fij and 
£j are provided by 



and 



-j E ( X ;-A;)( x /-Aj) 



with 9{j the number of vectors in the training set belonging to 
class CO;. The decision process then corresponds to assigning 
a vector x to class co 7 if gj (x) > g^ (x) for all k ^ j, where 

1 



S ; (x)=lnP(co 7 ) 



In 



(*-£;)> 



an estimate of P (co,-) being given by 



P(C0,) 



First we tested the classifier by using all the 1600 DFA vec- 
tors as the training set. This yields functions gk (x) that are 
able to correctly classify all vectors, a flawless performance. 
In order to evaluate the generalizability of the classifier re- 
sults, we randomly selected 80% (1280) of the 1600 available 
vectors to define the training set, using the remaining vec- 
tors in the testing stage, and took averages over 100 distinct 
choices of training and testing sets. When the training vectors 
were fed back to the classifier, as a first step toward validation, 
again no vectors were misclassified. Table |H| summarizes the 
average performance of the classifier when applied to the test- 
ing vectors. Notice that the maximum average classification 
error corresponds to 3%, for microstructure 8. The overall 
testing error, taking into account all classes, corresponds to 
0.8%. Misclassifications almost exclusively involve assign- 
ing a vector to a microstructure with the same average domain 
size but different average density (e.g. microstructures 4 or 
12 instead of 8), with very few cases involving the same den- 
sity but with a different although closest average domain size 
(microstructure 8 instead of 7). The low error rate is due to 
the fact that, although the average DFA vectors of similar mi- 
crostructures (such as 1 and 2) lie close to each other, the vari- 
ances of data along the different directions have a sufficiently 
distinct profile. 



Microstructure 


Success rate 


Misclassifications 


1 


98.6 (0.3) 


5: 1.4 


2 


98.7 (0.2) 


6: 1.3 


3 


99.6 (0.1) 


4: 0.1 7: 0.3 


4 


99.0 (0.2) 


8: 1.0 


5 


99.1 (0.2) 


1:0.4 9:0.5 


6 


99.8 (0.1) 


2:0.1 10:0.1 


7 


98.8 (0.2) 


3: 1.13 8: 0.04 11: 0.03 


8 


97.0 (0.4) 


4: 2.9 12: 0.1 


9 


99.5 (0.1) 


5: 0.5 


10 


99.8 (0.1) 


14: 0.2 


11 


99.6 (0.1) 


7:0.3 15:0.1 


12 


99.4 (0.2) 


8: 0.6 


13 


99.7(0.1) 


9: 0.3 


14 


100 


none 


15 


100 


none 


16 


99.1 (0.2) 


12: 0.9 



Table II: Average performance of the classifier when applied to the 
testing vectors, calculated over 100 sets of training and testing vec- 
tors. For the second column, bold numbers indicate the percentage of 
vectors which were correctly classified, figures in parenthesis corre- 
sponding to the standard deviations. The third column registers aver- 
age misclassifications (in the first row, '5: 1.4' indicates that 1.4% of 
vectors belonging to microstructure 1 were misclassified as belong- 
ing to microstructure 5). 



The efficiency of the classifier is dependent on the choice of 
values of the time-window sizes. For instance, restricting the 
values of T; to the powers of 2 doubles the overall testing error, 
to around 1.7%, while expanding } to the integer parts of 
the powers of 2 1 / 8 leads to a much larger overall testing error 
of 24%. Choosing |Ty} as the integer parts of the powers of 
2 1,/2 actually leads to a slightly smaller overall testing error of 
0.6%, but a few training errors also occur. Thus, our choice 
of {t^} from the integer parts of the powers of 2 1 / 4 seems 
to be close to optimal. In contrast, performing the detrended 
fluctuation analysis according to the original recipe of Ref. 1 8 1 
leads to a minimum overall testing error of 28%. 



V. NEURAL NETWORKS 



In the spirit of Refs. [22] and [23 1, which employed artifi- 
cial neural networks in order to identify disorder parameters 
in the random-bond random-field Ising model, we wish to in- 
vestigate if a similar approach can be useful in estimating av- 
erage domain sizes and average densities based on fluctuation 
analyses of our simulated ultrasound data. 

The idea here is to build a neural network which reads the 
DFA vectors as inputs, targets as outputs the physical param- 
eters (either average domain size or average density) from all 
vectors of 15 of the 16 possible microstructure, and then tries 
to guess the corresponding parameter from the DFA vectors 
of the remaining class. The network — a multilayer percep- 
tron 11241 l25l — is composed of an input layer with N\ = 37 
neurons, which receive the data from each DFA vector, an out- 
put layer with a single neuron (N4 = 1), whose reading is the 
desired physical parameter, and two hidden layers, contain- 
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Figure 5: Real (circles) and predicted (squares) values of the rescaled 
average domain sizes for different samples (numbered 1 to 1600). 




1000 

sample 

Figure 6: Real (circles) and predicted (squares) values of the rescaled 
average densities for different samples (numbered 1 to 1600). 



ing respectively N2 = 18 and N3 =8 neurons. The connection 
weights between neurons in contiguous layers are adjusted so 
as to minimize the mean square error between the desired and 
the actual outputs, according to the backpropagation prescrip- 
tion. We employed the hyperbolic tangent as the activation 
function, and both input and output data were converted to a 
logarithmic scale and adjusted so as to lie between —1+8 and 
1 — e, with e of order 1/10. (This rescaling improves the per- 
formance of the perceptron when dealing with microstructures 
for which parameters take extreme values.) 

Figures [5] and [6] show plots of the rescaled average domain 
size and rescaled densities for each microstructure, along with 
the predictions output when the perceptron is trained with data 
from all remaining microstructures. Despite the relatively 
high average error (of 0.05 for the domain sizes and 0.1 for 
the densities), it is clear that the information hidden in the 
DFA vectors is enough to provide useful predictions for the 
unknown parameters. We also trained the perceptron to out- 
put both domain size and density, but the performance showed 
a considerable decline compared to the case when the param- 




Figure 7: Sketch of the geometry used in simulating ultrasound prop- 
agation in an inhomogeneous medium. See main text for labels. 



eters where targeted by different networks. 



VI. SUMMARY 

Our aim in this work was to provide, within a fairly con- 
trolable framework, a proof-of-principle for the identification 
of microstructures based on fluctuation analyses of ultrasound 
signals. With a slightly modified detrended-fluctuation anal- 
ysis (DFA) algorithm, we were able to build an automated 
Gaussian classifier capable of assigning a DFA curve to the 
correct microstructure among 16 possibilities, corresponding 
to combinations of 4 average densities and 4 average domain 
sizes, with an error rate below 1%. Although not detailed 
here, an analogous classifier based on the original DFA al- 
gorithm of Ref. [8 1, despite not providing a comparable per- 
formance (yielding a much larger error rate around 30%), is 
able to separate the microstructures according to their average 
densities with an error rate of about 2%. Incidentally, yet an- 
other analogous classifier based on Hurst's R/S analysis also 
performs modestly for overall classification, with an error rate 
of about 22%, but is able to separate the microstructures ac- 
cording to their average domain size with a success rate in 
excess of 99.7%. 

We also described a multilayer perceptron which is able to 
provide estimates of a physical property for DFA curves from 
an unknown microstructure after being trained to output the 
corresponding property for the remaining microstructures. 

The application of the methods described here to more real- 
istic situations depends on a series of tests which incorporate 
effects coming from more complicated, higher-dimensional 
geometries. Among these, we mention mode conversion at 
domain interfaces and the presence of additional defects such 
as voids or inclusions of distinct phases. We hope the results 
reported in this paper will encourage future investigations. 
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Appendix A: Solution of the wave equation for a 
one-dimensional heterogeneous medium 

We want to solve the wave equation for the displacement 
field 3>(x,t), 



3 2 4> 

dx 2 



i a 2 * 



with the sound velocity c a constant along each domain into 
which the one-dimensional system, of size L, is divided. This 
means that <t> (x, t ) will be given by a different function in each 
domain, and the problem can be recast as the solution of the 
wave equations 



d 2 <t>j _ 1 d 2 <S>j 



dx 2 



c] dt 2 



for every domain j, subject to the boundary conditions 



and 



Pj-i c j-i 



2 3*y-i 



dx 



(xj,t) 



(Al) 



(A2) 



(A3) 



describing the continuity of the displacement and the pressure 
fields at the interfaces between domains. Here, Cj and pj de- 
note the sound velocity and the density in domain /', while xj 
is the coordinate of the left end of domain j. The medium 
is divided into N + 2 domains (j g { — 1,0, 1,2, . . . ,N}), with 
X-i =0 and xjv+i = £; see Figure [7] Domains from j = 1 to 
j = N correspond to the medium to be investigated, and span 
a length W < L. Domain j = holds a piezoelectric trans- 
ducer, in which the ultrasonic pulse is to be produced, and 
domain j = —1 is reserved for an "escape area", introduced 
to mimic the presence of an absorbing wall at the back of the 
transductor. 



Separation of variables leads to a general solution of Eq. 
(JaTJ, for a given angular frequency co, of the form 

<&j(x,t;(o) = \AjCos(kjx)+BjS\n(kjx)]cos((dt) 
+ [Cj cos (kjx) + Dj sin (kjx)] sin (cor) . 

Since we will impose initial conditions in which ^J- = 0, we 
can set Cj = Dj = for all /'. The boundary conditions in Eq. 
(|A"2]> and (fAl) lead to 



Aj-i cos (kj-iXj)+Bj-i sin (kj-iXj) 
Aj cos (kjXj) + Bj sin (kjXj) 



and 



z j- i\Aj-ism(kj-iXj)-Bj-i cos ixj)} 
Zj [Aj sin (kjXj) - Bj cos (kjXj)} 



= 1 



(A4) 



(A5) 



for j £ {0, 1,2,... ,Af} , in which we introduced the acoustic 
impedances Zj — PjCj, while the reflective boundary condi- 
tions at x — and x = L, 4>_i (0,f) = 4>at (L,t) = 0, yield 



(A6) 



A-i=0, 

Aa? cos (^L) + sin (^atL) = 0. 
In order to mimic an absorbing wall at the back of the trans- 
ductor (at x = xq), we choose for domain j = — 1 an acous- 
tic impedance z-\ = zo and a very small sound velocity c_i. 
These choices guarantee that waves incident on the left end 
of the transductor are not reflected, rather entering the escape 
area and not returning during the simulation. 

Equations ( |A4| i, ( | A5[ > and ( |A6| l constitute a homogeneous 
system of linear equations in the coefficients Aj and Bj. 
Rewriting the system as the matrix equation 



Aq Bo 



An B n )M w = 0, 



in which 



Ma 



1 cos(fe_ixi) z-i sin(fe-i^i) 

sin(fc_ixi) — z-i cos (k-\xi) 

- 

- 




cos (kox\) 
sin (koX\) 





-z sin (fen^i) 
zocos(^oxi) 





cos (kox\) 
sin (koX\ ) 

— cos (fclXl) 

— sin(A;ixi) 



zosin (^l) 
-zo cos (koxi ) 
-zisin 

Zicos(^ixi) 

























••• cos (kff—iXff) ZN-i$in(kN-iXN) 

••• sin (fey-l^Af) — zat_i cos (k^-iXfj) 

••• — cos (kfifXiy) —Zn sin (k^XN) cos (kpjL) 

••• — sin (JcnXn) zncos^nxn) sinlk^L) 



we see that nontrivial solutions for the |A ; -,Z?j} are obtained only if 



detMAf = 0, 
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an equation whose solutions correspond to the eigenfrequen- 
cies (ti k . Minor-expanding the determinant using the last col- 
umn of Mn leads to the recursion formulas 



with 



: detM„ = - cos (k n x n+ i) f„ + sin (k n x n+ \ ) g„, 



f„ = sin (k„x„)h n— 1 "F Zn COS (/^X^) 17l n — \ , 

g n = cos(k n x„)h n -i-z n sin(k n x n )m n -i, 
h„ = -z„sm(k n x n+ i) f n -z n cos (k„x„+i)g 



subject to the initial conditions 



m_i = sin(&_ixo) 



and 



h-\ = -z_icos(^_ixo). 



These formulas allow the numerical evaluation of the deter- 
minant for an arbitrary value of co. For a given geometry, the 
eigenfrequencies are numerically determined by first sweep- 
ing through the values of co, with a certain increment 8co, un- 
til detMw changes sign, bracketing an eigenfrequency whose 
value is then refined by the bisection method. The process 
is repeated with decreasing values of Sco until no additional 
eigenfrequencies up to a previously set value C0 max are found. 
For a given eigenfrequency co^, the corresponding coefficients 
{AjicBjjA (with a new index k to indicate the dependence on 
CO^) are determined recursively from Eqs. ( A4i and (A6 1, sup- 



plemented by the normalization condition = 1. 

The general solution of the full wave equation then takes 
the form 



4> (x, t ) = Y,§kXk (x) cos (()) k t ) . 

k 



in which 

X k (x) = A jk cos ((Okx/cj) + Bj k sin ((ti k x/cj) , 

with j such that Xj < x < Xj+\. The constant coefficients fa 
are derived from the initial condition <t> (x, 0) by using the or- 
thogonality condition satisfied by the X k (x), 

N f x i+i 

Y, Pj dxX k (x)X q (x)=^ k d kq . (A8) 

j=-i Jx i 

Explicitly, we have 



= r L Pi f dx<S> W . (A9) 

/=— 1 Jxj 



with ^ k > defined by Eq. (A8i. The conclusion that the 
orthogonality condition involves the densities p ; - comes from 
integrating, over the entire system, the differential equation 
satisfied by X k (x), multiplied by X q (x), then exchanging the 
roles of q and k, and subtracting the results, taking into acount 
the continuity of the pressure at the interfaces. The above 
orthogonality condition follows if K> k ^V) q . 

From the displacement field, the sound pressure can be cal- 
culated as 



p(x,t) = pjc 2 j 



d<&(x,t) 
dx 



■ PjcjY § kX k W C0S (®kt ) , 



(A7) again with j such that Xj <x < Xm. 
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